Melting and Rippling Phenomena in Two Dimensional Crystals with localized bonding 
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We calculate Root Mean Square (RMS) deviations from equilibrium for atoms in a two dimen- 
sional crystal with local (e.g. covalent) bonding between close neighbors. Large scale Monte Carlo 
calculations are in good agreement with analytical results obtained in the harmonic approximation. 
When motion is restricted to the plane, we find a slow (logarithmic) increase in fluctuations of the 
atoms about their equilibrium positions as the crystals are made larger and larger. We take into 
account fluctuations perpendicular to the lattice plane, manifest as undulating ripples, by examin- 
ing dual-layer systems with coupling between the layers to impart local rigidly (i.e. as in sheets of 
graphene made stiff by their finite thickness). Surprisingly, we find a rapid divergence with increas- 
ing system size in the vertical mean square deviations, independent of the strength of the interplanar 
coupling. We consider an attractive coupling to a flat substrate, finding that even a weak attraction 
significantly limits the amplitude and average wavelength of the ripples. We verify our results are 
generic by examining a variety of distinct geometries, obtaining the same phenomena in each case. 

PACS numbers: 62.25.Jk, 62.23.Kn,63.22.Np 
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I. INTRODUCTION 

Efforts to gain a quantitative microscopic understand- 
ing of melting have spanned more than a century. The 
Lindcmann criterion developed in 1910^ describes melt- 
ing in terms of the Root Mean Square (RMS) devia- 
tion from the atomic equilibrium positions. Since long- 
range positional order stems from the periodic arrange- 
ment of atoms in crystalline solid, atomic deviations that 
are comparable to the separation between atomic species 
could obscure the regularity of the underlying crystal lat- 
tice with a concomitant loss of positional order. The 
Lindemann criterion specifies that melting has occurred 
if the RMS deviations reach on the order of a tenth of 
a lattice constant, and has proved to be a reasonably 
effective theory for three dimensional systems. 

The Lindemann analysis does not take into account 
correlations of the motions of neighboring atoms. Corre- 
lations are more important at lower dimensions, and the 
process of melting is hence strongly dimensionally depen- 
dent. While three dimensional crystals exhibit long-range 
order below certain temperatures, statistical fluctuations 
play a significant role in one dimensional systems, pre- 
cluding all but short-ranged local ordering for T > 0. 

The process of melting in two dimensions is more sub- 
tle, and is understood in the modern context to occur in 
more than one stage. An initial continuous loss of posi- 
tional order precedes the proliferation of lattice defects, 
which accumulate and eventually complete the melting 
process at sufficiently high temperatures by destroying 
even orientational order, where each atom has a fixed 
number of neighbors. 

Thermally induced fluctuations in atomic positions can 
have an important effect on nano-engineered systems 
where structures may be on the atomic scale. Atomic 
clusters or quantum "dots" are mesoscopic assemblies of 
atoms where the scale is confined in all directions. Lin- 
ear structures such as carbon nanotubes are essentially 



one dimensional objects (although having cross sections 
on the atomic scale) where the tube length may ap- 
proach macroscopic scales. Finally, two dimensional sys- 
tems with nanoscale thickness such as covalently bonded 
graphene sheets are genuine monolayers with thicknesses 
on the atomic scale, but spanning macroscopic areas. 

The novel charge transport properties of graphene have 
been of intrinsic fundamental interest, and have also in- 
spired scenarios for the use of graphene in semiconduc- 
tor microprocessor applications. Technological uses for 
graphene will need a stable planar substrate for the im- 
plementation of nano-circuitry, and fundamental scien- 
tific research will also benefit from the minimization of 
the amplitude of random undulations in graphene layers. 

We examine two dimensional crystals with properties 
that would generically be found in two dimensional co- 
valently bonded crystals, including stiffness with respect 
to displacements perpendicular to the plane of the sheet. 
Although we do not consider temperature regimes capa- 
ble of disrupting the lattice topology or number of neigh- 
bors for each atom (e.g. by thermal rupture of bonds be- 
tween neighboring atomic species), we examine the loss 
of order caused by fluctuations of atomic positions about 
their equilibrium positions which nonetheless leave the 
bonding pattern intact. 

If the motion of particles comprising the crystal is con- 
fined to the plane of the lattice, the gradual loss of long- 
range crystalline order with increasing system size has 
been understood as being in some respects similar to the 
destruction of ferromagnetic ordering in the X-Y model 
(the motion of the spins are confined to the plane with a 
ferromagnetic coupling between them) by thermally ex- 
cited spin waves. Nevertheless, on a detailed level the two 
systems differ. In the case of the X-Y model, spin-spin 
correlation functions decay algebraically with the spatial 
separation between spins below the Kosterlitz-Thouless 
temperature for vortex unbinding. On the other hand, 
the RMS deviation in atomic positions in two dimen- 
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FIG. 1: Square lattice extended coupling geometry with in- 
teractions with nearest and next-nearest neighbors (left) and 
the periodic triangular lattice and labeling scheme (right) 



sional crystals has been described as logarithmically di- 
vergent (i.e. varying as log [a (T)L] where a(T) is a tem- 
perature dependent parameter) for any finite tempera- 
ture^. 

In Section II, we discuss theoretical techniques and the 
system geometries under consideration. Then, in section 
III we examine three dimensional lattices where we show 
directly for suitably rigid lattice geometries that the RMS 
deviations from equilibrium converge to a finite value 
in the thermodynamic limit, an anticipated property of 
three dimensional systems. Moreover, we determine a 
reference temperature threshold Tjf where mean square 
fluctuations about equilibrium reach one tenth of a lattice 
constant, corresponding to the melting point according to 
the Lindemann criterion. This temperature will serve as 
a point of reference in the examination of two dimensional 
systems where thermal fluctuations disrupt long-range 
order for any finite temperature. However, although we 
find stable crystalline order in three dimensional geome- 
tries, we also discuss a significant caveat which applies 
for simple cubic lattices and other geometries which lack 
local stiffness. To a great extent, the lattice geometries 
we report on are based on the two dimensional examples 
shown in Fig. [TJ A square lattice pattern, and a trian- 
gular lattice structure are shown. The former lacks in- 
herent rigidity, but the square lattice gains local rigidity 
through the activation of an extended coupling scheme 
in which both nearest neighbors and next-nearest neigh- 
bors interact. In the same way, a simple cubic lattice re- 
quires interactions between next-nearest as well as near- 
est neighbors to resist thermal fluctuations and maintain 
long-range crystalline order. 

By considering two geometries and appropriate three 
dimensional generalizations which differ in significant 
ways (i.e. one base on a square pattern and the other 
assembled of triangles or tetrahedra joined at their cor- 
ners), we identify generic thermodynamic characteristics 
common to both. 

In Section IV, we examine two dimensional lattices 
such as those shown in Fig. [1] with motion confined to 
the crystal plane, finding the very slow (i.e. logarithmic) 
loss of crystalline order anticipated for two dimensional 
crystals. On the other hand, for a two dimensional crys- 
tal embedded in three dimensions, it is important to con- 



sider transverse perturbations tending to push atoms out 
of the plane. We find that in the absence of binding to 
a substrate, two dimensional crystals are much less able 
to resist extraplanar distortions than fluctuations which 
are confined to the lattice plane. 

In Section V, we examine dual-layer systems where 
coupling between the crystal planes imparts local stabil- 
ity with respect to extra-planar variations of atomic po- 
sitions in a caricature of physical systems (e.g. sheets of 
graphene) which have a finite thickness, and would be im- 
bued with local stiffness. We find for two distinct locally 
rigid dual-layer geometries similar rapid divergences of 
mean square displacements as the crystal is made larger, 
corresponding to thermally induced rippling of the crys- 
tal, and scaling linearly with the size of the system. Anal- 
ysis of the density of vibrational states reveals that the 
length scale of the random undulations increases with the 
size of the system with strong long wavelength contribu- 
tions. On the other hand coupling to a flat substrate, 
however weak, places an asymptotic upper bound on the 
ripple amplitudes and also limits the average wavelength 
of thermally induced undulations. 



II. CALCULATION METHODS AND 
MONTE-CARLO SIMULATION RESULTS 

We examine thermodynamic properties (e.g. the mean 
square deviations of atoms about equilibrium positions) 
for crystals with short range bonding in the regime where 
bonds remain intact and thermally induced lengthening 
and shortening of bonds is small relative to the unper- 
turbed, or equilibrium, bond length. With individual 
bonds varying only slightly in length, it is appropriate 
to model the bonds as harmonic potentials so the cou- 
plings between neighboring atoms are effectively treated 
as springs connecting the two particles. It is important 
to note that although we neglect anharmonic effects from 
the bonds, the noncollinearity of bonds in the crystal ge- 
ometry may in principle introduce anharmonic terms in 
the Hamiltonian. Nevertheless, at temperatures near and 
below the melting point, many scenarios are amenable to 
the harmonic approximation where the neglect of anhar- 
monicities (whether intrinsic or geometric) has a small 
impact on the accuracy of the calculation. Analytical 
results obtained in the context of the harmonic approx- 
imation are validated in the cases we consider by good 
agreement with Monte-Carlo calculations where the an- 
harmonic characteristics of the bonding stemming from 
peculiarities of the lattice geometry are rigorously taken 
into account. 

The Hamiltonian is given by 
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where is the equilibrium energetically favored bond 
length and Z,j is the instantaneous separation between 
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atoms i and j. The outer sum is over the atoms in the 
(finite) crystal, and the inner sum is over the neighbors 
associated with the atom indexed by the label i. The 
additional factor of 1/2 is included to compensate for 
double counting of bonds. The constant is the sec- 
ond derivative of the interatomic potential Vy (r) at the 
equilibrium separation . 

We develop the harmonic approxima- 
tion directly from the bond length Zy = 
yj(xi - xj) 2 + (yi - yj) 2 + (zi - Zj) 2 For the x co- 
ordinates, it is convenient to write, for example, 
Xi = x® + Sf where xf is the equilibrium coordinate and 
Sf is the shift about equilibrium. We operate in the 
same way for the y and z coordinates, finding 



(A°f + 5f-s*) 2 + (^y + sy-5«) 2 



where A?f 
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,,.-„ _(»?-!,»), and A?/ = 
z®). One may develop the harmonic approximation by 
expanding terms such as (Zy — Zy) to quadratic order in 
the shift differences (Sf-Sf), (Sf-Sf), and (5? -6$). The 

result will be (Zy — If, ) w Ay ■(£,• — £,) , where Ay 



is a unit vector formed from Ay = (A??, A°J, A? z ).The 
terms <5i and 5j are vector atomic displacements such 
that, e.g., Si = (Sf, S?, Sf). A salient characteristic of the 
bond energy is its dependence on the differences of the 
coordinate shifts (e.g. Sf — SJ for the x direction) instead 
of Sf , S\ , and Sf by themselves, a condition which under 
many circumstances permits the neglect of anharmonic- 
ities due to bond non-collinearity. 

In the harmonic approximation, the lattice energy due 
to deviations from equilibrium positions will be 
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On expanding, one obtains a quadratic expression mixing 
the displacements 
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Diagonalizing the appropriate matrix yields 3iV eigen- 
vectors, taken to be normalized. Each of the set of 3N 
eigenvectors has a component for the individual degrees 
of freedom in the crystal lattice, permitting the lattice 
Hamiltonian to be written in decoupled form as 



K 3N 

T E AqC ' 



(5) 



with eigenvector expansion coefficients c a and eigenval- 
ues A a ; the parameter K is the "primary" harmonic con- 
stant, which is taken to be the nearest neighbor intra- 
planar coupling constant in schemes, such as extended 



models with multiple coupling constants. The eigen- 
modes are independently excited by thermal fluctuations, 
and thermodynamic equilibrium observables may be cal- 
culated by evaluating Gaussian integrals. As an exam- 
ple, the thermally averaged mean square fluctuation per 
atomic species (<5rms) is (first moments of the coordinate 
shifts such as (Sf ) vanish in the thermal average and do 
not appear in the expression below) 

N 

<<W} 2 = £ E^) 2 + (S\f + (Sf) 2 ) (6) 
i=i 

Indexing the eigenvectors with the label a and noting, 

3N 

e.g., that Sf = ^^c Q t>", we see that the total square of 

a=l 

the instantaneous fluctuations per particle is 

N 3N 3N 

w = ^EE E +«) 

(7) 
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In calculating the thermal average the term c Q c Q ' will be 
as often negative as positive when a/a, and there will 
only be a non-zero contribution to (<5r,ms) 2 ii a = a ■ 
Hence, the double sum enclosed in square brackets will 
collapse to a single sum, and the calculation is reduced 
to the thermal average 

3N N 

<<w 2 = £ E< c2 >E [^) 2 + (^) 2 + «) 2 ] ( g ) 

a—l i—l 

The eigenvector normalization condition gives 



E[(^) 2 + (<) 2 + (^) 2 ] = 1 > 
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(3) and hence (<5rms) 2 appears simply as 
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The partition function Z may be calculated with the aid 
of jT ) 00 e ~ aq dq = (tt/o) 1 / 2 , and one has a product of 
decoupled Gaussian integrals, which may be written as 
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with j3 = 1/fcs, &b the Boltzmann constant, and the 
temperature T is given in Kelvins. For the sake of con- 
venience, units are chosen such that the lattice constant 
a is equal to unity, and a reduced temperature is defined 
with t = kftT/K . Evaluating the integrals in the product 
given in Eq. [IT] yields for Z 
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The thermally averages mean square displacement may 
be written in terms of a thermal logarithmic derivative 
of Z, and in particular, one finds 
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Hence, the thermally averaged mean square deviation 
from equilibrium may be written as the square root of 
a sum over eigenvalue reciprocals. 
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Zero eigenvalues would lead to a diverging expression, 
but eigenvalues which are strictly equal to zero are ar- 
tifacts of periodic boundary conditions, correspond to 
global translations of the crystal lattice, and are excluded 
from the sum. The dependence on reduced temperature 
consists of a i 1 / 2 factor. To concentrate on characteristics 
specific to a lattice geometry and its coupling scheme, as 
well as trends with respect to system size L, the nor- 
malized mean square displacement <52 MS will often be 
discussed in lieu of the full temperature dependent quan- 
tity. 

In the case of a periodic regular crystal lattice, it is 
useful exploit translational invariance, which will lead to 
exact expressions for the vibrational mode eigenstates 
and frequencies for periodic crystals (or at the very least 
yielding a small matrix which may be diagonalized an- 
alytically or by numerical means if necessary) if atomic 
displacements are written in terms of the corresponding 
Fourier components. 

Using Monte Carlo calculations to sample thermody- 
namic quantities incorporates anharmonic effects in a rig- 
orous manner, providing a means of assessing the validity 
of the harmonic approximation. We employ the Metropo- 
lis technique 3 - to introduce random displacements and 
sample the distribution corresponding to thermal equi- 
librium. We follow the standard Metropolis prescription, 
where an attempted random displacement with an as- 
sociated energy shift AE is accepted with probability 
e -AB/fc B T if AE > q and the Monte Carlo 

move is in- 
variably accepted for cases in which AE < 0. 

In calculating thermodynamic quantities, we operate 
in terms of Monte Carlo sweeps where a sweep, on av- 
erage, consists of an attempt to move each atom in the 
crystal with the acceptance of the move subject to the 
Metropolis condition. In the calculations, the sampling 
of thermodynamic quantities is postponed until the com- 
pletion of the first 25% of the total number of sweeps to 
eliminate bias from the initial conditions, which are not 
typical thermal equilibrium configurations for the sys- 
tem. To reduce errors due to statistical fluctuations in 
the Monte Carlo simulation and obtain several digits of 
accuracy in the results, we conduct at least 5 x 10 5 sweeps. 
Figure[5](for the square lattice with an extended coupling 
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FIG. 2: Graph of the RMS fluctuation about equilibrium 
versus systems size L for various values of the reduced tem- 
perature t for the square lattice with next-nearest neighbor 
couplings. The solid lines are analytic results obtained in the 
harmonic approximation, and symbols are results from Monte 
Carlo calculations. 



scheme) and Figure [3] (for the triangular lattice geome- 
try) show mean square deviation curves for various tem- 
peratures ranging from temperatures an order of magni- 
tude smaller than to temperatures on par with the 
Lindemann criterion result for the melting temperature 
of the corresponding three dimensional system. The solid 
lines correspond to analytical results, while the symbols 
are RMS values obtained with Monte Carlo calculations. 

The curves show very good agreement between the 
Monte Carlo data and analytical results over a wide range 
of temperatures and system sizes, and deviations are pri- 
marily mild statistical errors (on the order of one part in 
10 3 ) in the Monte Carlo calculations. 



III. RIGID AND NON-RIGID THREE 
DIMENSIONAL LATTICES 

To establish a temperature scale for the two dimen- 
sional systems, where long-range crystalline order is not 
expected to exist at temperatures above OK, we first ex- 
amine three dimensional lattices, which may exhibit long- 
range positional order at finite temperature if the lattice 
is suitably rigid. As a preliminary step, we perform an 
analysis similar to the Lindemann treatment where an 
atom in a simple cubic geometry is coupled to six nearest 
neighbors. Since we do not take into account the motion 
of neighboring atoms, we take their displacements to be 
zero; certainly the excursions of neighboring atoms would 
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FIG. 3: Graph of the RMS fluctuation about equilibrium ver- 
sus systems size L for various values of the reduced temper- 
ature t for the triangular lattice. The solid lines are analytic 
results obtained in the harmonic approximation, and symbols 
are results from Monte Carlo calculations. 



Prom the Gaussian integration, we find Z = (7rf) 3 / 2 . The 
RMS displacement will be (r 2 ) 1 / 2 where 



t 2 jrLn{Z) = 3t/2 



(18) 



Hence, the thermally averaged mean square shift is 
(St/2) 1 / 2 . The Lindemann criterion places the melting 
temperature at a temperature high enough that the mean 
square deviation <5rms reaches a tenth of a lattice con- 
stant, which corresponds to tjP = |(10~ 2 ), a reduced 
temperature on the order of 0.01. 

If correlations among atoms are taken into account, 
next-nearest neighbor couplings become crucial to im- 
parting local stiffness and maintaining long-range crys- 
talline order. To see how rigidity is an important factor, 
we calculate the RMS displacements for a simple cubic 
lattice where only couplings between nearest neighbors 
are taken into account. The energy stored in the lattice 
will be 



K n_1 

e = - y 

2 ^ 

i,j,k=0 



-(4 



ijfc+1 



- 8 V 

ij-\-lk ijk 



(19) 





Simple cubic geometry 
(not rigid) 



Next nearest neighbor 
coupling included 
(rigid configuration) 



FIG. 4: Illustration of the simple cubic, nonrigid structure 
and rigidity gained by incorporating next nearest-neighbor 
couplings as shown in the image to the right. 



average to zero, although to be precise, one would need 
to take into account cooperative effects of the atomic mo- 
tions of the neighbors. 

The lattice energy has the form 
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where, for example, A x+ is the shift in length of the 
bond to the nearest neighbor in the positive x direc- 
tion. Applying the harmonic approximation and taking 
the atomic shifts to be {6 X , 6 y , S z }, the energy becomes 



K[5l 



(16) 



In the calculation of 5rms j the partition function has the 
form 



Z 



— oo J — oo J — oo 



' dS x dS y dS z 



(17) 



where a periodic geometry is assumed, and the counting 
factor of 1/2 does not appear since the sum has been 
constructed to avoid redundancies. If we use the trans- 
formations 



fix e I(k x i+k y j+k z k) 



(20) 



where I is the imaginary unit, and similar expressions 
are used for 8j- k and 5^ k . In terms of the Fourier com- 
ponents, the lattice energy may be written as 



(1 — cos k 
+ (1 
+ (1 



■)K\ 2 
)K\ 2 



cos k 7 



(21) 



The x, y, and z degrees of freedom 6%. , S^, and 5^. auto- 
matically decouple. The normalized mean square devi- 



ation is yj ~^2 a A Q 1 , where the sum is restricted to non- 
zero eigenvalues. We identify three eigenvalues, a£ = 

2(1— cosfcx), = 2(1 — cosky), and A^ =2(1 — cosfc z ) 
for each wave vector {k x , k y ,k z } As can be seen in Figure, 
the mean square fluctuation about equilibrium positions 
grows very rapidly with increasing system size. The di- 
vergence in the RMS displacements is a consequence of 
the lack of rigidity of the simple cubic geometry, which 
facilitates the destruction of long range crystalline order 
by thermal fluctuations. However, next-nearest neighbor 
couplings make the lattice rigid, and are very effective 
in suppressing fluctuations about equilibrium and estab- 
lishing long-range crystalline order for the simple cubic 
lattice. 

The structure of the eigenvalue density states profile 
has informative characteristics particular to the lattice 
geometry from which it is obtained, and the density of 
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FIG. 5: RMS displacements graphed with system size L 
for various values of K2, expressed in units of the nearest- 
neighbor coupling Ki. Panel (a) shows a closer view of the 
5rms curves over a smaller range of system sizes, and panel 
(b) is a graph with a broader range of system sizes included 
in the plot. 
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FIG. 6: RMS displacements graphed as a function of the ratio 
of next nearest neighbor to nearest neighbor coupling, K% /Ki . 
Inset (a) is a graph of the normalized RMS deviation with 
respect to log 10 (-K2/-Ki), and inset (b) shows the RMS fluc- 
tuations raised to the 1/4 power versus log 10 (_K2 / Ki) ■ 
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FIG. 7: Normalized Eigenvalue Density of States for the 
simple cubic system for an extended coupling scheme with 
K 2 = {0.0,0.01,0.02} with a sampling of 2.4 x 10 s eigenval- 
ues. 
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states is calculated for many of the systems we report 
on. We achieve the thermodynamic in a genuine sense 
by not restricting k x , k y , and k z to discrete values as is 
done for finite systems. The density of states is built up 
by Monte Carlo sampling in which the wave- vector com- 
ponents are each generated independently from a uni- 
form random distribution. To obtain good statistics, at 
least on the order of 2 x 10 s eigenvalues are sampled in 
constructing the DOS. The same Monte Carlo sampling 
procedure is used to calculate the <5rms values shown in 
Fig. |51 and thereby completely remove any bias from fi- 
nite size effects. The density of states corresponding to 
the simple cubic system (shown in the graph in Fig. [7]) is 
consistent with the divergence of the RMS fluctuations 
with increasing system size. The bimodal structure is 
sharply peaked in the low and high eigenvalue regimes, 



FIG. 8: Normalized Eigenvalue Density of States for the 
simple cubic system for an extended coupling scheme with 
K 2 = {0.05, 0.10, 0.20} with a sampling of 2.4 x 10 s eigenval- 
ues. 



with the former contributing to the steady rise of <5rms 
with increasing system size L. 



In the extended coupling scheme in the simple cubic 
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FIG. 9: Normalized Eigenvalue Density of States for the 
simple cubic system with coupling for an extended coupling 
scheme with K 2 = {0.5, 0.7, 1.0} with a sampling of 2.4 x 10 8 
eigenvalues. 



geometry, the energy stored in the lattice is 
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(22) 



where K\ is the coupling to nearest neighbors, K 2 is the 
coupling to next-nearest neighbors, and k 2 = K 2 /K\ is 
the ratio of the next-nearest and nearest neighbor cou- 
pling constants. In terms of Fourier components, one has 
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with k indicating the wave- vector with components k Xl 
k y , and k z , and again k 2 — K±/K 2 . The eigenvalues are 



hence obtained by diagonalizing the 3x3 matrix 
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«2 sin fc z sin fc^ k 2 sin fe 2 sin fc y 
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(24) 



Although solving the cubic characteristic equation will 
yield analytical expressions for the eigenvalues, the re- 
sult is cumbersome, and we instead use standard algo- 
rithms for the diagonalization of a symmetric matrix to 
efficiently obtain the eigenvalues numerically. 

The eigenvalues determined in this manner are used to 
calculate the means square atomic fluctuations, and the 
results are shown in Figure [SJ where <5rms is graphed with 
respected to L for a range of the next to nearest neigh- 
bor coupling strength ratio k 2 — K 2 jK\. Whereas the 
mean square displacement steadily rises with system size 
when k 2 = (i.e with only nearest-neighbor couplings ac- 
tive), the curves behave very differently for nonzero k 2 , 
ultimately saturating with increasing L. The stabiliza- 
tion of <5r,ms in the thermodynamic limit indicates the 
presence of intact long-range crystalline order. In the 
case of k 2 =0, the mean square deviation steadily di- 
verges with increasing L. The same divergence with only 
nearest-neighbor interactions taken into a account occurs 
whether one is considering the simple cubic structure, a 
square lattice, or a linear chain. Hence, in some lattice 
geometries, having a three dimensional structure may be 
insufficient to stabilize long-range order if an extended 
coupling scheme is not taken into consideration. 

By switching on and varying the strength of the next- 
nearest coupling K 2 , one sees the appearance of long- 
range crystalline order as the cubic system is made in- 
creasingly rigid. In Fig. [5] the mean square displacement 
is shown graphed versus the coupling ratio K 2 /K\. The 
tendency for atoms to be driven from their positions in 
the lattice does increase as K 2 is shut off, but the diver- 
gence occurs at a slow rate. Inset (a) is a graph of <5rms 
versus the logarithm of K 2 jK\. While the concavity of 
the curve indicates a somewhat more rapid than loga- 
rithmic divergence, a semi-logarithmic plot of <5rms (i- e - 
as shown in inset (b) of Fig. [5]) shows an approximately 

linear scaling of £ R ms with the logarithm of the system 
size, still a relatively slow divergence, albeit somewhat 
more rapid than a simple logarithmic divergence. Hence, 
the next-nearest neighbor couplings in the extended cou- 
pling simple cubic model are very effective in restoring 
long-range crystalline order. 

Trends in the eigenvalue density of states profile with 
increasing K 2 / K\ to next-nearest neighbors are shown in 
in Fig. [71 Fig. [5J Fig. [5] The almost immediate retreat 
of the low and high frequency peaks toward the center is 
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consistent with the effectiveness of an extended coupling 
scheme in stabilizing long-range crystalline order even 
for very small values of the ratio ft 2 — K%jK\. The DOS 
profile has a simple structure for small K2 < 0.1, while 
intermediate K2 values are associated with a richer den- 
sity of states curve which changes rapidly as the coupling 
ratio is increased further. 

As in the case of the cubic lattice with the extended 
coupling scheme, we may calculate the lattice energy in 
the harmonic approximation for the tetrahedral lattice, 
which is 



x-{Si + ijk—5ijk 

(±x+^y)-(4+ifc-4^ 
^■y)'(Si+ij-ik-Sijk) 



k ™ _1 

e=- y 

i,j,k=0 



fa 

I- 



lz)-(Sijk+i-5ijk) 



2V3 



y+ 



(Si- 



ljk+l—0ijk. 



y+viz) • (Sij-ik+i—Sijk) 



(25) 



where there is only one coupling constant K since bonds 
are considered between nearest neighbors only, the tedra- 
hedral geometry being intrinsically rigid, and we have 

used the fact that the altitude of a tetrahedron is 



times the lattice constant. The energy may be expressed 
in terms of Fourier components, and one has the task of 



FIG. 11: Normalized root mean square (RMS) deviation 
shown versus log 10 L for the three dimensional tetrahedral 
crystal. The inset is a graph of the normalized RMS devi- 
ations, again plotted with respect to log 10 L, with the hori- 
zontal line indicating the extrapolated <5rms in the thermo- 
dynamic limit. 



diagonalizing the 3x3 matrix 
/ 4 - 2cosfc a; - \ 



I (cos ky+COS k z ) 
2 cos(/by k x ) 
\-\cos,(k z -k y )J 



^fcos(k y -k x ) 
2 \ — cos fcj 



2 (cos(k z — k. 
3\ — cosk z 



^fcos{k y -k x 
2 V — cos k„ 



( 4 



"2 cos ky 




"2 cos(^ky kj<^ 
— ^ cos k z — 

cos(k z -k z )J 



2cos(k z — k y ) 
— cos k z — 




2icos(k z - k x 
3 \ — cos k 7 



(2 cos(k z — k y f\ 
- cos k z - 
cos(k z — k 2 



3 — cos k z — 



cos(k z — k x ) I \cos(k z — k y ) 



(26) 

The three dimensional tetrahedral lattice is locally stiff 
even with only nearest-neighbor couplings taken into ac- 
count, and the rigidity inherent in the tetrahedral lat- 
tice geometry is sufficient to preserve long-range crys- 
talline order, as may be seen in Figure [TT] which displays 
the normalized mean square deviation versus the system 
size L. The inset is a semi-logarithmic plot with the 
horizontal axis extending over three decades of system 
sizes. The saturation of the normalized RMS displace- 
ment with increasing L is evident in both of the graphs, 
and in the thermodynamic limit, is in the vicinity of 1.12. 
With temperature dependence included, one will have 
<5rms = 1.12t 1 / 2 . Hence, the Lindemann criterion would 
give tg D = 0.0080, compatible with the previous estimate 
which neglected correlations of the atomic displacements 
from equilibrium. 
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yielding the eigenvalues 
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FIG. 12: Normalized Eigenvalue Density of States for the cu- 
bic crystal with nearest and next-nearest neighbor couplings 

(a) , and the density of states profile for the tetrahedral lattice 

(b) . 



In inset (a) of Fig. Q21 the density of states is shown 
for the simple cubic lattice geometry with the extended 
coupling scheme, and for the tetrahedral lattice in inset 
(b). For both systems, while other details of the density 
of states profiles differ, the curves tend swiftly to zero 
in the small eigenvalue regime, a hallmark of intact long 
range crystalline order in rigid three dimensional lattices. 



IV. INTRAPLANAR MOTION 

We first examine the case where motion perpendicu- 
lar to the plan is forbidden, and atomic deviations from 
equilibrium are confined to the lattice plane. We consider 
various geometries, but first we consider a square (ef- 
fectively a face-centered system) , illustrated in Figure [U 
where coupling to the four next nearest neighbors is taken 
into account. We then consider a locally rigid triangular 
lattice where each atom interacts with six nearest neigh- 
bors. In both the face-centered square and triangular 
systems, we find a logarithmic divergence with increas- 
ing system size in the mean square fluctuations about 
equilibrium. 

For the periodic square geometry with the coupling 
scheme extended to next-nearest neighbors, the lattice 
energy to quadratic order is 



E = -Y 
2 ^ 



[,(, 



i+lj ' 



-S x - 



+ 



i+lj+l " 



(27) 



Operating in reciprocal space, one diagonalizes the 2x2 
matrix 



2 — cos k x \ . , . , 

, , sin k r sin ft,, 
cos k x cos kyl y 



. , . , / 2 — cos k v 

sin k T sin ft,, , " , 

y cos k x cos k y 



4 — cos k x — cos k y — 2 cos k x cos k y ) (29) 



±y (cos k x — cos k y ) + 4 sin k x sin k y 

In the case of the triangular lattice with six fold coor- 
dination, one may also obtain analytical expressions for 
the mean square deviations. In real space, the harmonic 
approximation for the energy stored in the lattice is 



K 

2 ^ 

hj=0 



x- 1 S i+ ij - Sij 



V?, 



Sij+l — Sij 



\ 



J 



(30) 



Expressing the displacements in terms of Fourier compo- 
nents, one decouples the x and y degrees of freedom by 
diagonalizing the matrix 



3 — 2 cos k x - 

i COS ky — 7; cos [ft 



V 



— 2 — (cOS^Aly 



- COS fc„ 



^ (cosffty — k x \— cosky) 



yields the eigenvalues 



3 q cos ky 
~2 cos[Aiy k x 



(31) 



= [3 — cosfta; — cosky — cos(k y — k x )} (32) 

/ COS 2 fca; + cos 2 /c y + COS 2 ( fc j, — k x )~ cosk x cos k y 
y —cosky cos(k y — k x ) — cos(k y — k x ) cosk x 

For convenience in comparison with the analytical re- 
sults in the harmonic approximation, we consider peri- 
odic boundary conditions in the crystal plane. We have 
also examined anchored lattices, where atoms at the pe- 
riphery are prevented from moving, while those in the 
interior are free to move. For both the free and fixed 
boundary conditions, as in the three dimensional case, we 
obtain qualitatively similar results, and the same physi- 
cal phenomena. 

In Fig, [inland Fig.[TH the normalized mean square de- 
viation <5^ M g is graphed with respect to the system size L 
for the square lattice in the extended scheme and the tri- 
angular lattice, respectively. The overall behavior of the 
mean square deviations from equilibrium is qualitatively 
the same for both lattice geometries. In both cases, the 
main graph is semi-logarithmic with (<5r MS ) 2 on the ordi- 
nate. The traces are linear to a very good approximation 
for all regimes of L (i.e. for small, moderate, and large) 
shown, and the linearity is maintained for four decades 
of system sizes ranging from several to on the order of a 
few times 10 4 lattice constants. 
(28) In Fig- E2 an d Fig. [TH inset (a) is a standard plot, 

and the apparent saturation of the <5^ M g curve is a hall- 
mark of the slow loss of long-range crystalline order best 



10 




0.5 



1.5 2.5 3.5 

Log 10 ( L) 



4.5 



FIG. 13: Square of the normalized root mean square (RMS) 
deviation shown versus log 10 L for the square lattice system 
with extended couplings. The solid line encompassing the 
open circular symbols is a strictly linear fit. Inset (b) is a 
semi- logarithmic graph of the normalized RMS deviations, 
plotted with respect to log 10 L. 



seen on a semi-logarithmic graph. Inset (b) in Fig. [13] 
and Fig. [5] contains as semi-logarithmic plot with S^ MS 
[ instead of (^ Mg ) 2 ] on the ordinate axis. The curves 
plotted in this manner are not linear, and it is evident 
that the divergence of the fluctuations about equilibrium 
is actually somewhat slower than logarithmic; instead, it 
is (<5r M s) 2 which scales as Ln(L). 

As in the case of the three dimensional systems, it 
is informative to examine the density of states, shown 
in the graph of Fig. [15] for the square lattice in the ex- 
tended coupling scheme in panel (a) and the triangular 
lattice in panel (b) of Fig. [15] Again, while details of 
the density of states profiles shown are peculiar to the 
lattice under consideration, the behavior in the regime 
of low eigenvalues is quite similar, and both curves tend 
to a finite value instead of dropping swiftly to zero as in 
the density of states for the rigid three dimensional lat- 
tices. The failure of the density of states to vanish in the 
small eigenvalue limit contributes to the slow divergence 
of <5rms in L. 



V. EXTRA-PLANAR MOTION 

The locally stiff face-centered square and triangular 
lattices show the anticipated slow logarithmic divergence 
in system size. However since laboratory systems often 
are not vertically constrained, it is important to exam- 
ine a scenario where motion perpendicular to the plane 
of the lattice may be considered. There is an impor- 
tant difficulty with single layer systems, in that motion 
perpendicular to the plane is not hindered since there 
are no restraining bonds with a directional component 




4000 8000 12000 
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Log 10 (L) 

FIG. 14: Square of the normalized root mean square (RMS) 
deviation shown versus log 10 L for the triangular lattice. The 
solid line encompassing the open circular symbols is a strictly 
linear fit. Inset (a) is a standard plot of the RMS devia- 
tion with respect to system size L, while inset (b) is a semi- 
logarithmic graph of the normalized mean square deviations, 
plotted with respect to log 10 L. 
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FIG. 15: Normalized Eigenvalue Density of States for the face 
centered square lattice with motion confined to the lattice 
plane, depicted in panel (a), and for the triangular lattice in 
panel (b). 



transverse to the plane of the layer. 

However, by considering dual-layer geometries, it is 
possible to incorporate local stiffness with respect to per- 
turbations that would push atoms above or below the 
lattice. We examine analogs of the simple cubic lattice, 
where we again use an extended coupling scheme to cre- 
ate local stiffness. On the other hand, we also consider 
a dual-layer tetrahedral lattice. Although the two lat- 
tice geometries achieve local stiffness in different ways, 
the similarities we find in thermodynamic behavior of 
the mean square atomic fluctuations suggest these char- 
acteristics would appear in the generic case as well. 

Fig.[l6]illustrates the structure of the dual- layer square 
lattice with an extended coupling scheme; the additional 
couplings between next nearest neighbors impart local 
stiffness to the system with respect to perturbations per- 
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FIG. 16: Illustration of the periodic dual-layer square lattice 
with nearest and next-nearest neighbor coupling and labeling 
scheme; blue indices refer to the upper layer, while red indices 
pertain to the lower plane. 




FIG. 17: Schematic representation of a coarse-grained super- 
structure for a graphene sheet 



pendicular to the planes of the square lattices. 

Fig-HHshows how the dual-layer systems is constructed 
as a caricature of the graphene lattice. The image la- 
beled (a) is a schematic illustration of a single hexagonal 
cell in a graphene monolayer. The bonding shares sim- 
ilarities with that in a benzene ring with delocalized tt 
orbitals forming honeycomb networks of charge density 
above and below the plane occupied by the carbon atomic 
nuclei. The superimposed lattice work is a rigid network 
compatible with the symmetries of the graphene layer 
and set up to capture the rigidity of the hexagonal cells 
making up a sheet of graphene. With the honeycomb 
graphene pattern removed, the remaining lattice geom- 
etry and the labeling scheme for the crystal members is 
shown in Fig. [T5] 

With the superscript I representing the lower plane 
and II indicating the upper plane, the lattice energy in 



FIG. 18: Illustration of the periodic dual-layer triangular lat- 
tice and labeling scheme; blue indices refer to the upper layer, 
while red indices pertain to the lower plane. 



the harmonic approximation is 



2 ^ 

i,j=G 



E 



+ K Z < 



A0+*)-(%H-%)] a 
3s(-*+*)-(*5-i-%) 



(33) 



The corresponding complex Hermitian 6x6 matrix to 
be diagonalized has the form 



aixix ai x iy aixiz ai x \i x a\ x \\ y a\ x u z 

dlylx a lyly aiyl Z aiyUx aiyUy O^TIz 

aizix a\ z iy ai z i z a\ z i\ x aiziiy aiziiz 

O-Uxlx O^Uxly O-llxlz QllMIx a l\xl\y llMb 

dllylx a llyly CLllylz a UyIIx CLUylly OlIj/IIz 

Mizix a\iziy auziz aiiziix aiiziiy auziiz 



A B 
fit A 



(34) 



where A and B are 3x3 matrices and W is the Hermitian 
conjugate of B. The sub-matrices A and B are given by 



A =K 



K z +4-2cosfc a A . . 

- , , 2 sm k T sm k„ 
—2 cos k x cos kyj y 







2 sin k x sin k y 




; z + 4 — 2 cos k % 
-2 cos fc T cos fc, 



3k, 



(35) 
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where k z = K z /K for A and 
cos k x 



B =K 7 



-i sin fc T 







cos k,. 



-isin fe 



isink x — isinfcj, — (l+cos/cz + cosfcj,). 



(36) 



for B where the six eigenvalues for each wave-number pair 
are calculated numerically with code available in the EIS- 
PACK linear algebraic library for diagonalizing complex 
Hcrmitian matrices. 

On the other hand a locally stiff dual-layer system 
which may be regarded as a section of the three dimen- 
sional tetrahedral lattice such that the upper and lower 
layers are triangular lattices with connections between 
the layers. The dual-layer lattice structure based on the 
tetrahedral geometry is illustrated in Figure IT51 The ver- 
tices of the upper layer are positioned above the centers 
of the triangles in the lower layer with bonds extending 
from atoms in the upper layer to each of the corners of 
the triangle below such that each atom in the dual-layer 
system is a member of a rigid tetrahedron; the result 
is a locally stiff layer, as in the dual-layer square lat- 
tice extended model, but with a very different geometric 
structure. 

The lattice energy for the dual-layer tetrahedral system 

c -> 1 r - 11 

to quadratic order in the displacements Sij and Sij has 
the form 



E 



K 



' i,j=0 



[Hk+ij-stif- 



\ 



[(§-^)-fe-i-%) 



+ 



y 

2s/3 

4=4- 



§*V0§-% + i) 



+ 



(37) 



Expressing the lattice energy in terms of Fourier com- 
ponents leads to a 6 x 6 matrix to be diagonalized, which 
may again be written in terms of 3 x 3 submatrices as 

A B 

~. -> , where 



f + 3 — 2 cos k x \ 
~ ^2 cos(/ty k% 
— ^ cos ky 



^3 / cos(k y ~k x ) 
2 I — cos k y 



^/3 / cos(k y -k x ) 



2 I 



cos k„ 



K * h3 

3 
2 



3 

cos (A: y 



2 ..' " 2 cos (^y) | () 







2k 2 



(38) 



where again k z = K z /K for the sub-matrix A and 



_( 1+e -ife«) _^ e -ik x _q v /|( e - l fe x _ 1 ) 



4 



V3 ie lj 3l+3e- lfc ' 



3 I e -tfc B 



1- 



-.—ik 



—e 



(39) 



for the sub-matrix i3. 

The use of a dual-layer lattice geometry to provide re- 
sistance to transverse deviations is not sufficient to pre- 
vent a rapid divergence in <5rms with increasing system 
size L. Whereas thermally averaged mean square fluctu- 
ations grew very slowly (i.e. logarithmically) when the 
atomic motions are confined to the lattice plane, <5rms 
for the dual-layer systems increases linearly with L; ulti- 
mately, it is not difficult for the sheet to bend and flex in 
the presence of thermal fluctuations in spite of its locally 
stiff characteristics. The diverging mean square devia- 
tions from equilibrium and other thermodynamic char- 
acteristics of the the dual-layer square lattice in the ex- 
tended coupling scheme and its counterpart based on a 
tetrahedral geometry are examined, with consideration 
given to the effects of increasing L and variations in the 
inter-layer coupling strength. 

The graph of the normalized RMS displacement, 
shown in Fig. [T5]for the dual-layer square lattice with an 
extended coupling pattern shows a dependence on sys- 
tems size which is an asymptotically linear growth in L. 
8rms curves for several values of the interlayer coupling 
K z are shown; the intra-layer coupling is taken to be 
unity, so K z is effectively expressed in units of K . Al- 
though the relative interlayer couplings range over three 
orders magnitude, there is little variation of the curves, 
especially for K z = {0.1,1.0,10.0}. Similarly, the <5rms 
curves ultimately vary linearly in the system size L with 
little dependence on the relative magnitude of K Zl which 
again is expressed in units of K . 

Again, it is useful to examine the density of states for 
the eigenvalues in the case of the locally rigid dual-layer 
systems, which are richer than the density of states pro- 
files corresponding to rigid three dimensional lattices or 
those of the single layer geometries with atomic fluctua- 
tions confined to intra-planar motion. Although details 
of the density of states profiles for the two geometries 
differ, the both curves show a divergence of the density 
of states with decreasing eigenvalue magnitude, whereas 
the density states remained constant in the case of the 
planar systems with exclusively intra-planar motion and 
vanished altogether for the rigid three dimensional sys- 
tems. Inset (a) of Fig. [2T1 show the DOS for the dual-layer 
square lattice, while inset (b) is a graph of the density of 
states for the locally stiff tetrahedral system. The DOS 
cusp for both lattices at the zero eigenvalue point is re- 
sponsible for the rapid divergence of <5rms with increasing 
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FIG. 21: Normalized Eigenvalue Density of States for the 
dual-layer cubic system with with nearest and next-nearest 
neighbor interactions (a) and for the dual-layer locally stiff 
lattice based on a tetrahedral lattice geometry (b) for system 
size L — 5001. 



FIG. 19: Normalized mean square displacements for K z — 
0.01, 0.03, 0.1, 1.0, and 10.0 for the dual-layer square lattice 
with an extended coupling scheme, where K z is in units of 
the intralayer coupling K. The inset is a closer view of the 
RMS curves. 
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FIG. 22: Normalized density of states profiles for K z — 3.0, 
1.0, 0.3, and 0.1 for the dual-layer square lattice with an ex- 
tended coupling scheme where K z is in units of the triangular 
lattice intra-layer coupling K. 



FIG. 20: Normalized mean square displacements for K z — 
0.01, 0.1, 1.0, and 10.0 for the dual-layer tetrahedral lattice 
where K z is in units of the triangular lattice intra-layer cou- 
pling K. The inset is a closer view of the RMS curves. 



L. 

While adjusting the interaction between the layers to 
enhance the resistance to local transverse perturbations 
has little effect on the mean square fluctuations for large 
values of L, the eigenvalue density of states evolves as 
the interplanar to intraplanar coupling ratio k z = K z /K 
is modified. Density of states profiles for k z values rang- 
ing from k z = 0.1 to k z — 3.0 are shown for the dual- 
layer square system with and extended coupling pattern 
in Fig. [22] and for the tetrahedral counterpart in Fig. [23J 

Density of states profiles are shown for strong (k z = 



3.0) and moderate (k z == 1.0) values of the the coupling 
ratio in insets (a) and (b) of Fig. [22] and Fig. [231 and 
there is little change in the DOS curve in the low eigen- 
value regime. On the other had, as k z decreases further 
and the interplanar coupling begins to fall below parity 
with that in the plane, the eigenvalue density of states 
profiles begin to change more drastically, as may be seen 
in panels (b) and (c) of Fig. [22] for the dual-layer square 
system and Fig. [23] for the dual-layer tetrahedrally based 
geometry; the distribution in both cases rapidly grows 
narrower with decreasing k z . Although the two lattice 
geometries are very distinct, similar (and likely generic 
to locally rigid dual-layer lattices) trends may be seen in 
the DOS profiles in the regime of low eigenvalues as n z 
is reduced. 
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FIG. 23: Normalized density of states profiles for K z — 3.0, 
1.0, 0.3, and 0.1 for the dual-layer tetrahedral lattice where 
K z is in units of the triangular lattice intra-layer coupling K. 



VI. COUPLING TO A SUBSTRATE 




System Size, L 



FIG. 24: For the dual-layer square lattice in the extended 
scheme, the main figure and the inset are graphs of mean 
square fluctuations versus system size L for a range of sub- 
strate couplings K 3 , with K a given in units of the inter-planar 
and intra- planar coupling constant (both equal to K). The 
symbol legend on the main plot also pertains to the inset 



We incorporate an attractive interaction with a flat 
substrate by including an additional harmonic poten- 
tial acting on the lower members of the tetrahedral and 
square extended coupling dual-layer systems. We take 
the attraction to depend only on the atomic shift <5?J 
above the planar system, and the additional term hence 
has the form fy. (^) 2 

Figure [24] shows the effect of the substrate coupling 
on Sums m the case of the dual-layer square system with 
an extended coupling scheme. The graph, which shows 
mean square deviation curves for a wide range of K a 
values, indicates the capacity of even a very mild sub- 
strate coupling to suppress thermally induced undula- 
tions in the dual-layer sheet. Similarly, for the tetra- 
hedrally based dual-layer lattice geometry, an attractive 
interaction with a substrate considerably reduces fluctu- 
ations transverse to the lattice planes, preventing a rapid 
divergence of £rms- The mean square deviation curves 
are shown in Figure 1251 

We also examine the effect of an attractive substrate 
coupling on the density of states profiles, and results are 
displayed in Fig. [26] for a range of substrate coupling 
constants K s . With increasing K s , a salient trend is the 
opening of a separation between the sharp cusp and the 
zero eigenvalue mark on the abscissa. The migration of 
the maximum formerly at the zero eigenvalue point to 
a peak at a larger eigenvalue is associated with a sharp 
reduction of the mean square fluctuations about equilib- 
rium, and the lattice is better able to withstand trans- 
verse fluctuations. 

The presence of a flat substrate plays a very impor- 
tant role in dictating the overall structure and ampli- 
tude of ripples in the dual-layer geometries we report on 



here. This result is in accord with recent experiments on 
graphene sheets deposited on cleaved mica substrates 4 , 
where the careful preparation of flat substrates signif- 
icantly dampens the ripple amplitude, whereas much 
larger undulations are seen with sheets attached to sub- 
strates with poorer control over the flatness. 

To determine which length scales are associated with 
the strongest contributions to the thermally averaged 
mean square deviations about equilibrium, we have pre- 
pared histograms showing the relative contribution to 
<5rms versus inverse wave- vector magnitude, with the lat- 
ter providing a length scale. Apart from a significant 
diminution in the height of thermally excited undula- 
tions in the dual-layered sheets, we also find a consider- 
able reduction in their typical wavelength. In figure 1271 
for the dual-layer tetrahedral lattice in the absence of a 
substrate couplings, the dominant contribution to i5rms 
comes from large length scales comparable to the scale of 
the lattice. However the picture changes with the activa- 
tion of a finite substrate coupling as may be seen in the 
inset with the peak height at the minimal wave-number 
decreasing with increasing K s . Moreover, as may be seen 
in Figure 1281 introducing even a weak anchoring to the 
foundation below immediately creates a strong peak in 
the short wavelength regime, skewing the size of ther- 
mally induced ripples toward smaller length scales. 

In conclusion, we have examined thermally induced 
fluctuations about equilibrium in two and three dimen- 
sional crystalline solids with a local bonding scheme. 
While long-range crystalline order may exist in three 
dimensional crystal lattices, some geometries (e.g. the 
simple cubic lattice) are not rigid when only nearest- 
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FIG. 25: For the dual-layer tetrahedral lattice, the main figure 
and the inset are graphs of mean square fluctuations versus 
system size L for a range of substrate couplings K a , with K s 
given in units of the inter-planar and intra-planar coupling 
constant (both equal to K). The symbol legend on the main 
plot also pertains to the inset 



| 0.8 
£ 0.7 

0.6 
£0.5 
J 0.4 
•a 0.3 

1 0.2 
g 0.1 
I 0.0 



(a) 




Ke =0.0 





0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 

Eigenvalue Magnitude 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 

Eigenvalue Magnitude 




1.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 

Eigenvalue Magnitude 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 

Eigenvalue Magnitude 



FIG. 26: Normalized Density of States for the Dual layer 
system for various substrate coupling strengths K s . Panels 
(a), (b), (c), and (d) show the density of states for K s equal 
to 0.0, 0.1, 0.5, and 1.0 respectively. 



neighbor couplings are taken into account, and an ex- 
tended coupling scheme is needed to prevent the diver- 
gence of mean square fluctuations with increasing system 
size L. 

In two dimensional lattices, we find RMS fluctuations 
to increase at a very slow (logarithmic) rate when motion 
is confined to the lattice plane. On the other hand, when 
transverse motion is permitted, thermal fluctuations are 
very effective in bringing about significant vertical dis- 
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FIG. 27: Normalized Eigenvalue Density of States for the 
dual-layer triangular system with system size L = 5001. 
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FIG. 28: Normalized Eigenvalue Density of States for the 
dual-layer triangular system with system size L = 5001. 



placements of particles which contribute to rapidly grow- 
ing deviations from equilibrium, and 5rms ultimately 
diverges at a linear rate in L. The asymptotically lin- 
ear divergence in the mean square deviations from equi- 
librium is insensitive to the strength of the interlayer 
coupling; <5rms values appear to converge and eventu- 
ally show identical behavior with increasing system size 
whether the coupling K z established between the layers 
to provide local stiffness is quite weak or very strong rel- 
ative to the bonding K between atoms in the same layer. 

Introducing a coupling K s to a flat substrate very effec- 
tively hinders transverse fluctuations in two dimensional 
crystal lattices, even in the coupling is very weak, and re- 
flects the importance of a substrate in shaping the char- 
acteristics of ripples set up by thermal fluctuations by 
inhibiting transverse deviations. An attractive coupling 
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to a fixed substrate also reduces the typical lateral length induced ripples, 
scale or wavelength of thermally excited undulations in 
lattices bound to a substrate. These tendencies are con- 
sistent with recent experimental observations that con- 
trol over the flatness of the underlying surface is directly Acknowledgments 
related to the amplitude and length scale of thermally 
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